Eigenfrequencies of the randomly pinned drum and conductivity of graphene 
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Graphene is convenient material for nanomechanichal applications since high-frequency oscilla- 
tions are easily accessible. In this Article, we consider graphene on a rough substrate attached to 
imperfections at random locations. We explore the statistics of low-lying phonon modes, which 
exert most influence on the conductivity of graphene. We find that the nearest neighbor spacmgs 
of low lying eigenfrequencies have the Wigner-Dyson probability distribution after averaging over 
the random configurations of disorder. Due to interaction of electrons with the oscillations of the 
membrane, an electron can be transfered to higher or lower energies, which is a manifestation of 
the phonon-assisted Tien-Gordon effect. The Tien-Gordon effect suppresses the conductivity of 
graphene. In the regime of low Fermi energies and small sizes of the sample an increase of con- 
ductivity is observed which we refer to Klein tunneling and electron-hole pair creation. Eventually, 
when the increase of the transmission becomes too prominent, the pair creation changes the ground 
state of the system, signalizing the limit of applicability of the single-particle Dirac equation used 
in this paper. 



I. INTRODUCTION 

Graphene is a monolayer of graphite with excel- 
lent elastic properties and low mass, which make it a 
prospective nanomechanical resonator— attractive for a 
wide range of applications from mass sensing in single- 
molecule ranged to quantum manipulation of the elemen- 
tary mechanical vibrations^. They are combined with 
perfect electrical properties. The mobility of the carri- 
ers, both in suspended samples^ and the samples on the 
substrate^, is one of the highest among all modern semi- 
conductors. The combination of mechanical and electri- 
cal properties opens wide opportunities for manipulation 
and read-out of mechanical oscillations using electrical 
signals. 

In most experiments, graphene is deposited on sub- 
strate. The substrate is often not ideal and can be con- 
sidered as a random potential landscape. The graphene 
membrane interacts with the substrate, which may re- 
sult to its attachment to the substrate at some areas or 
points, and free suspension between these areas. This ar- 
rangement was confirmed experimentally^. This means 
that the mechanical oscillation spectrum of such mem- 
brane is complex, with localized and extended phonons 
determined by the substrate profile. Moreover, since the 
profile is random, the frequencies of these excitations can 
be considered as random quantities. It is known that the 
frequencies do not uniquely characterize the shape of the 
oscillating membrane. Indeed, this question was stud- 
ied in 1966 by Marc Kac in his famous paper "Can one 
hear the shape of the drum"—, with the conclusion that 
is possible to construct different in shape, but isospectral 
membranes. Therefore the statistical approach remains 
the only meaningful way to characterize the frequency 
spectrum of this system. 



In this Article, we characterize statistical properties of 
these frequencies and show how they can be assessed via 
electric transport measurements. 

The elastic energy of the oscillating membrane has two 
contributions: stretching and bending energies. We con- 
clude that the bending energy is much smaller than the 
stretching one. Hence the dynamics of graphene mem- 
brane is described by the Hclmholtz's equation with the 
fixed boundary conditions at the regions where the mem- 
brane is attached to the substrate. The same equations 
describes the oscillation of the classical membrane. In 
our work we show how to " hear" the disorder of attach- 
ments of the membrane and see the distribution of the 
low frequencies averaged over the disorder. The frequen- 
cies can be directly measured in the experiment, but for 
graphene it is more convenient to extract information on 
the oscillation of the membrane from the conductivity 
measurements. 

The high-energy spectrum of the eigenvalues of the 
Laplace operator is well studied and is known to obey 
the Wigner-Dyson distribution^. This statement was 
proven by a variety of methods, analytical (random ma- 
trix theory and non-linear sigma-model) as well as nu- 
merical (various models of quantum billiards) . Therefore 
in this Article we study the probability distribution of 
the low energy spectrum of the membrane on disordered 
substrate. This question was not previously addressed 
since most of the studies of eigenvalues concern electron 
systems, where the relevant energy levels lie around the 
Fermi surface and have very high energy. In contrast, 
the phonon frequencies have a clear cut-off at low en- 
ergies. This cut-off frequency determines the scale for 
low-frequency modes. The second reason why low en- 
ergy levels are beyond the scope of the usual methods 
is that no statistical analysis can be performed for only 
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few levels. Thus, for a clean system the problem of the 
statistics of low-lying eigenvalues can not be formulated. 
One needs a disordered system to address this issue: Av- 
eraging over the disorder acts as averaging over the sta- 
tistical ensemble. The problem can be addressed only 
numerically as the analytical methods mentioned above 
are not applicable for low energy spectrum. According to 
our numerical modeling the distribution of low frequency 
eigenvalues of the membrane randomly attached to the 
substrate also reveals Wigner-Dyson statistics. 

Electrons in graphene are coupled to the deformation 
of the lattice. Due to the deformation of the lattice the 
overlap between the electron orbitals is changed, so that 
the band structure which describes the effective motion of 
the electrons is altered. Therefore the electron conductiv- 
ity can be influenced by the oscillations of the membrane. 
In monolayer graphene, the deformations can be strong, 
which leads to strong electron-phonon coupling. The cou- 
pling of the out of plane phonons to the electrons results 
in electrostatic potential and pscudomagnctic field which 
depend on the deformation, acting on electrons^. Since 
electrons move in the oscillating potential created by the 
phonons, they can scatter with absoprtion or emission of 
energy from the membrane. If the current-voltage char- 
acteristic of the system has some structure at the voltage 
V, then this structure propagates also to the voltages 
V + nftiUJi, with LUi being the set of the membrane fre- 
quencies and n an integer number. This phenomenon is 
similar to the Tien- Gordon effect which arises if the one 
lead of the sample is irradiated by the photons^. In the 
situation we consider, however, the oscillating potential 
is created by the membrane itself, and the area where the 
electrons feel this potential is separated from the leads by 
potential barriers. Then the conductance exhibits Fabry- 
Perot resonances at the energies Ei. Due to the oscilla- 
tions with frequencies uji each peak is enhanced at the 
energies Ei + nfuoi. If the conduction modes of graphene 
at different frequencies do not interact, the presence of 
the addition to the Fabry- Perot peaks provides us with an 
opportunity to measure the distribution of eigenfrequen- 
cies of the oscillating membrane. This way of detecting 
the oscillations of the membrane will work for narrow 
graphene strips. For wide strips, the mixture of Fabry- 
Perot peaks from different harmonics and their satellites 
will not give a clear picture. 

A number of papers^— studied phonons in graphene. 
They were focused on phonon spectrum either not influ- 
enced by disorder or influenced by weak disorder. We 
consider the special case, when the disorder changes 
macroscopic motion of the membrane, namely, the 
graphene sheet is attached to the substrate at random 
points. 

This Article is structured as follows. In Section HT1 we 
present the detailed results on the statistics of low-lying 
eigenfrequencies of the graphene membrane. The numer- 
ical method used to obtain the results is found in Ap- 
pendix [5J In Section IIIII the conductivity modification 
due to the Tien- Gordon effect is considered first for one 



frequency, then for two frequencies, and finally the in- 
fluence of the distribution of the phonon modes on the 
conductivity is addressed. In Section [IV] we discuss the 
results. In Appcndix|B]we show the method used to solve 
the scattering problem for determining the conductivity 
of the oscillating membrane. 

II. PHONONS IN DISORDERED GRAPHENE 

Graphene on corrugated substrate can be modeled as 
a membrane which is randomly attached to some regions 
of the substrate. We describe this membrane using the 
boundary integral method and confirm the usual wisdom 
that the distribution of the spacings between the near- 
est neighbour phonon energy levels for high levels is the 
Wigner-Dyson distribution. Moreover, we investigate the 
statistics for the low energies and find that, after aver- 
aging over disorder, it is also described by the Wigner- 
Dyson type function. 

A. Model 




FIG. 1: Left: Schematic representation of graphene attached 
at the arbitrary regions to the substrate. Graphene can os- 
cillate between the regions of attachment. Right: The model 
used to calculate the frequencies of oscillations of graphene 
membrane: circular membrane with circular regions of at- 
tachments placed randomly. 

Scanning tunneling miscoscopy (STM) experiments 
show^ that graphene on a rough substrate is attached 
to the substrate at some regions, at the hills of the sub- 
strate, and it is freely suspended over the valleys. To 
take this into account, we model a graphene sheet on a 
substrate as a membrane which is fixed at some regions. 
Due to external forcing or due to interaction with envi- 
ronment the membrane can vibrate. In this Section, we 
are interested in the statistics of eigenfrequencies of such 
membrane. 

We assume that the membrane is attached to the sub- 
strate strongly enough, so that it can not be detached due 
to oscillations. The depinning of the membrane on the 
rough substrate was considered in Ref. [IH The authors 
concluded that for the large scale of the imperfections of 
the substrate the membrane indeed remains attached to 
the substrate, namely for h/h < h^/S c / K(7/h) 2 , where 
h is the height of the fluctuations, I is the characteristic 
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length of the fluctuations, at which the attachment oc- 
curs, S c is the coupling strength of the graphene to the 
substrate, estimated for graphene in SiC>2 as ~ 2meV— , 
k is the bending rigidity of the two-dimensional layer, 
K ~ leV, and the length s cale h is the effective thickness 
of graphene layer h = n/ E2D, h « lA, E2D is the 
two-dimensional Young modulus, — 340A/m. For 
example, if the height fluctuations are of the order of 5 A, 
the attachment region of lnm is sufficient to keep the 
membrane attached to the substrate. We assume below 
that the condition is met, and depinning does not occur. 
The Lagrangian of the oscillating membrane is 



Pa4> 2 



(1) 



with cj)(t) being out-of-plane deviation. The first term 
in Eq. (JTJ stands for the kinetic energy, pg is the 
density of the two-dimensional graphene layer, po = 
7.6 x 10 -7 kg /m 2 , the second term is the bending en- 
ergy. The third term is the energy of the deformed mem- 
brane with the tension S, which can be estimated via the 
Hookc's law S = E2dAL/L with the relative elongation 
AL/L = 5o/E 2 d +£/L. Here S is the initial pretension 
and £ is the deformation of the layer (the amplitude of 
the oscillations). Hence we get the ratio of the bending 
and strectching energies E bend /E str ~ (h/L) 2 (AL/i) _1 . 
For a large membrane, this ratio is small, since, h/L <C 1. 
In the following, in order to obtain the eigenfrequencies 
of the membrane, we only take into account the term 
determined by stretching. 

Assuming that the oscillations are harmonic, cj> oc 
exp(iwi), we obtain the equation for the spectrum of the 
membrane with the border Q (at which the out-of-planc 
displacement vanishes) , 



A<f>(q) 



A 2 <j>{q) = 0, for q in n/dfl , 
4>(q) = for q in <9S1 , 



(2) 
(3) 



with A 2 = pqlu^/5. The boundary fi consists of the ex- 
ternal border of the graphene sheet fio and of the ar- 
eas {Hi, Q2, ■ ■ ■ , Qn} where graphene is attached to the 
substrate. We note that formally A 2 correspond to the 
energy in the equivalent Schrodinger equation widely dis- 
cussed in the context of quantum chaos. Therefore when 
we consider the statistics of the eigenfrequencies of the 
oscillations we always refer to A 2 , not to A„. 

The eigenvalues of the system @, are found nu- 
merically by using the boundary integral method-^, see 
Appendix [A] for more detail. This method is based on the 
re-writing the eigenvalue problem for the Laplace equa- 
tion © and its boundary condition in the form of 
the integral equation for the normal derivative of <f> on 
the boundary. It leads to solving the eigenvalue problem 
for the system of linear equations in a numerical imple- 
mentation of the method. We generalized the method 
for the case of several boundaries (see Appendix |A"|) . We 
tested the program on exactly solvable problems such as 
an eigenfrequencies of a circle and of a circle with a hole 
in the middle. 



To model graphene on substrate we take a circle with 
circular regions of the attachments in the middle (to be 
referred below as impurities or holes). The model is gen- 
eral enough as its classical analogue is non-integrable al- 
ready if two circles are randomly place inside the large 
circle. The impurities are placed at random positions. 
The eigenfrequencies are computed for each configura- 
tion. The averaging over disorder configurations is per- 
formed to determine the statistics. 



B. Statistics of eigenfrequencies at high energies 
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FIG. 2: The distribution of the nearest neighbour levels spac- 
ings is shown for the integrable case of the circular membrane 
and for the case of the membrane with holes. The distribu- 
tions are taken for the same range of the energy levels. For 
the integrable circular membrane the distribution (A) reveals 
the Poissonian distribution function @, N = 444. We note 
that the distribution is quite noisy. There is no ensemble av- 
erage for the integrable system, hence for the fixed range of 
energies there is a fixed number of levels. In the inset we show 
the distribution for larger range of energy, N = 7180, which is 
well fit by the Poissobian distribution. For the chaotic mem- 
brane (B) the distribution of the high energy levels reveals 
Wigner-Dyson probability distribution J5). The distribution 
is computed for the two randomly placed holes. The number 
of the levels is N = 401442. 
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Without any regions of the attachment in the middle 
the system (a circle) is integrable, and its eigenfrequen- 
cies obey Poissonian statistics, 

P(x) ~ exp(-z), x = X 2 n+1 - X 2 n . (4) 

see Fig. [2)3. A circle with two impurities in the middle 
becomes classically chaotic, namely a small initial dif- 
ference between trajectories grows infinitely in time. It 
is well-known that the statistics of high-energy levels of 
classically chaotic systems (such as non-integrable bil- 
liards) is well fit by the same distribution function^^. 
The same statistics appears in the disordered solid state 
problem treated by non-linear sigma-model, see Ref. [l9l . 
Indeed, we find that for the high levels the statistics of 
the nearest level spacings is Wigner-Dyson, 

P( s ) = | S exp(-f S 2 ) , (5) 

with s being the normalized spacing between two levels, 
see Fig. H\. 

To get higher eigenvalues, we need to discretize all 
boundaries very accurately, and it becomes computation- 
ally consuming to implement more holes. This is why we 
have chosen the minimal number of impurities (two) for 
which the system is chaotic. 

We note that if one models the attachment regions 
as points, the system is classically integrable, hence the 
statistics of the levels for a correctly chosen window of the 
energies is Poissonian^, even though for a smaller win- 
dow it is Wigner-Dyson. In this manuscript, we remain 
in the Wigner-Dyson regime. 

C. Statistics of low eigenfrequencies of the 
membrane 

The statistics of high-energy levels for our system is not 
different from any other system described by the Laplace 
equation. What is particular in a membrane over sub- 
strate is that we now have an access to statistics of low- 
lying levels as well, whereas for electron levels, which typ- 
ically are the subject of investigation in the level statis- 
tics problems, low- lying levels are not accessible and are 
outside the scope of the problem. For high levels, the 
random matrix theory appears as collective dynamics of 
the levels, and the statistics is effectively described by 
the Wigner-Dyson distribution. For low levels we do not 
have such averaging over the level positions. To the best 
of our knowledge, statistics of low-lying levels has not 
been studied so far since the absence of large number of 
levels means that there is no ensemble over which one 
might perform averaging to reproduce the results of ran- 
dom matrix theory (which is where the Wigner-Dyson 
distribution comes froro^i). It is not obvious what dis- 
tribution the lowest lying levels obey. 

We study this statistics in more detail than the one 
of high energy levels for different numbers of impurities. 
We treat the impurities as quenched disorder, hence we 




FIG. 3: The density of eigenfrequencies for the first few os- 
cillation modes of the membrane. The density is smeared by 
disorder in the position of the holes, but the oscillating behav- 
ior reflecting the positions of the eigenvalues in the absence 
of the disorder is present. We notice that in the region of the 
fourth-fifth initial eigenvalue of the membrane the density of 
states is already substantially more smooth that for the first 
few levels. The data is shown for two different sizes of the 
holes, r = 0.1 and r = 0.05. The averaging is done for the 
different number of the holes N = 5, 6, 7. 



always have average over disorder. The density of states 
averaged over different disorder configurations oscillates 
revealing the position of zeros of the clean membrane 
and becomes constant for large energies, see Fig. Ill CI 
We averaged over configurations with different numbers 
of the impurities, TV = 5, 6, 7. The oscillation in the den- 
sity of levels are more pronounced for smaller size of the 
impurities (the number of impurities is the same). It 
is in agreement with the simple-minded argument that 
small impurities have less influence on the spectrum. We 
argue that that the oscillations in the density of low en- 
ergy states do not disappear for the case of non-intcgrablc 
boundary as the oscillations simply reveal the geometric 
shape of the membrane. Therefore for a non-integrable 
outer boundary of arbitrary shape the oscillations would 
be shifted according to the first few energy levels of the 
membrane without disorder. 

We obtain the Wigner-Dyson distribution for the level 
spacings of the first four-five low lying energy levels 
slightly depending on the number of impurities in the 
system, Fig. We average all obtained numbers for 
different number of impurities and get the distribution 
which is close to the Wigner-Dyson one for different sizes 
of the impurities, Fig. [4j\,B. Let us note that we fit the 
numerical data with the dependence which is the same as 
in Wigner Dyson probability distribution. In the Wigner- 
Dyson probability distribution for high energy levels, the 
parameter s of the Wigner-Dyson probability function 
is normalized by the average level spacings <5A 2 , namely, 
s = (A 2 l+1 — A 2 )/5X 2 . In the case of low levels, the mean 
level spacing for the region of fit oscillates, see Fig. Ill CI 
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FIG. 4: Statistics of the low-lying eigenfrequencies. A: The 
distribution of the NNL (next nearest levels) for different 
numbers of impurities, N = 5, 6, 7. The distributions has sim- 
ilar properties as we see from the histograms and from the fit. 
B: The distribution of the NNL for the first few eigenfrequen- 
cies for different size of the impuritites, r = 0.1 and r = 0.05 
(averaging over different number of impurities, N — 5,6,7). 
We notice that it reproduces well the Wigner-Dyson distribu- 
tion function (J5). The normalization is done for the averaged 
to the constant density of states. 



the fit in Fig. @] is close to the Wigner-Dyson fit with the 
parameter s normalized by the mean value of density in 
the region of fitting. 

Furthermore, we explore the statistics of the position 
of the first energy level in the chaotic disordered cav- 
ity, Fig.EJ 

First we compare the distribution of the energy levels 
for one impurity obtained by solving the corresponding 
equation numerically and computed from the Born ap- 
proximation. We show the difference between conven- 
tional way of computing the shift of eigenenergies due 
to the disorder for the electron system with an impurity 
of finite potential V and the phonon system with fixed 
to zero boundary condition on the circumference of the 
impurity. 

In the Born approximation, the shift of the energy of 
the first level for the impurity Sl^ cdplaced at the random 
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FIG. 5: The energy distribution of the first eigenfrequency, 
Pjv(Ai), for different numbers of the attachement regions, N. 
A: Pi(\f) in comparison with the distribution in the Born 
approximation ©, (0. B: Evolution of the energy distribu- 
tion increasing number of holes, Pi,2,7(Ai). C: The energy 
distribution becomes stable for the large number of holes, 
JWCA?). 



position at the distance p from the center is 

E 1 (p) E[ Q) = V [ J 2 (^ 0) V)^ ? : « Vnr 2 jZ(E[ 0) p) , 

(6) 

where V is the potential characterizing the degree of at- 
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tachment of the membrane to the substrate, and Jo is 
the Bessel function of the zeroth order, which represents 
the solution of the Laplace equation in a circle with the 
zero boundary conditions. The probability to find the 
impurity at the distance p from the center is 



phonon motion of the membrane ai 



,11 



P{p) = Airrp. 



(7) 



We can get the dependence p{E\) from Eq. © and con- 
sequently the probability distribution P\ (E\ ) in the Born 
approximation. Fig. [5K demonstrates that the Born ap- 
proximation does not correctly describe the problem with 
the zero boundary conditions. The discrepancy can be 
attributed to the fact that we need an infinite potential 
to reproduce zero boundary conditions at the finite cir- 
cumference of the attachement region, while the Born 
approximation is only applicable for Vr 2 -C 1. For this 
reason, we proceed with numerical solution for a more 
complex system, with two and more impurities. 

As we would expect from simple-minded considerations 
analogous to those in the quantum mechanics, the eigen- 
frequency can move only up in energy with increasing 
the amount of disorder (amount of impurities). Thus, 
the minimum eigenfrequency corresponds to the one of 
the ideal circle. On all our plots we display the distri- 
bution of the first energy level taking into account this 
finite energy offset. The distribution of the first level 
is smoothened with increasing the disorder. In Fig. [5)3 
we show the evolution of the distribution with increasing 
number the impurities from N = 1 to N — 2 and N = 7. 
We clearly see how the distribution becomes more broad. 
In Fig. [5]C it is seen that the distribution becomes more 
or less similar for larger number of impurities, N = 5, 6, 7 

Microwave experiments^ are an established way to 
measure the eigenfrequencies for the wave equation, sys- 
tem given by Eqs. |(5J), ([3]), in various geometries. There- 
fore we can propose studying microwave radiation in the 
cavity as a direct way of measuring the computed dis- 
tributions. The boundary condition is imposed by plac- 
ing highly reflective material. One can think of a simple 
model: a big cavity with randomly placed disks. While 
changing the positions of the disks and measuring the 
first few eigenfrequencies, the distribution for the aver- 
age density and the spacing between nearest neighbors 
can be recovered. In this system, the lowest-lying eigen- 
value statistics must be the same at the one we described 
above. 



III. PHONON-ASSISTED TIEN-GORDON 
EFFECT 

In this Section, we investigate how the low-lying eigen- 
frequencies of the oscillating graphene membrane can be 
experimentally assessed by means of conductance mea- 
surement. 

The electrons traversing graphene are coupled to the 



Vl u _ ( 9l{UxX + Uyy) g 2 f*{Uij) \ 

V ^ h -\ g 2 f( Uij ) gi (u xx +u yy ) ) ' W 

where g\ is the deformation potential, g\ ~ 20 — ZtSeV , af- 
fected by screening; 172 is the coupling to pseudomagnetic 

2U X y + i(U XX Uyy). FOT 



fields, t/2 ~ 1.5 eV, and /(iiy) 



the second-order process involving out-of-plane phonons 
one obtains 



V. 



el—ph 



9i 92 



9i 
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9l c 



(9) 



where the similarity sign shows that we do not take 
into account the angle of the scattering between different 
phonon modes. We can make the following estimations of 
the parameters in Eq. q n ~ n/L, g\ ~ lOeV (includ- 
ing some screening). The amplitude of the oscillations <fi 
can be estimated in several different ways: as the ampli- 
tude of zero-point motion, <f) = (S/(2Mw)) 1 ' 2 , or as the 
amplitude of thermal fluctuations, Mcu 2 (/) 2 = KgT, or it 
can be varied by changing the amplitude of the oscilla- 
tions to which one pumps the energy (for example, by 
placing the sample on the piezosubstrate). 



Generalized Landauer formula for Tien-Gordon 
effect 



Reflected p Transmitted 

E + 2hw 




E + hu 



E 



E-hu . 



E - 



{<*} 
















— 








Vlead 
















Tien-Gordon effect 

FIG. 6: Left: While propagating through the region with po- 
tential oscillating with the frequency cj an electron with en- 
ergy E can scatter to the energies E + nhoj. Right: Schematic 
representation of the current flowing from left to the right 
through the graphene membrane oscillating with the fre- 
quencies {oji}. The Tien-Gordon effect happens due to the 
electron-phonon coupling. 



Let us consider electrons in graphene subject to oscil- 
lating potential. Electrons can emit or absorb quanta 
of energy corresponding to the oscillations due to the 
electron-phonon coupling. Hence an electron with the 
energy E which is coming to the region with a harmoni- 
cally oscillating potential with frequency u can be trans- 
mitted and reflected to the energies E ± nHui. A similar 
effect of the transmission of the electron from the energy 
E to the energy E + nfko due to microwave irradiation 
of the sample is called the Tien-Gordon effect. Note 
however that there is an important difference, since in 
the Tien-Gordon effect the bias voltage oscillates. Our 
system is formally equivalent to the conductor with the 
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oscillating barrier between the leads. Electrons in the 
leads are in thermal equilibrium described by the distri- 
bution functions /l{E) and fn(E + eVb), with Vb being 
the difference in the potentials between two leads. 

The electron propagating from the left to the right 
can be transmitted from the energy E\ to the energy 
E 2 or reflected to the energy E 2 - We denote the corre- 
sponding transmission and reflection coefficients from as 
Tei^e 2 an d Rei->e 2 - F° r the electron propagating from 
the right to the left the similar coefficients are denoted as 
Te 1 ^e 2 an d Re!^e 2 ■ Transmission and reflection proba- 
bilities are determined from the solution of the scattering 
problem similarly to Ref. I23I (see Appendix IB")). The total 
transmission from the left to the right is determined by 
the following expression: 

£ f dE (f L (E) R eU^eME + ifiu,)- 

where q is the transverse momentum. The first term rep- 
resents electrons moving to the right. The second and the 
third terms represent reflected electrons moving back to 
the left and electrons transmitted from the right, respec- 
tively. 

The change of integration variable gives: 

/ dE RE+thu^Eh{E + ihw) = 
J i 

= I ' dE^RE-tE-itohiE). (10) 

i 

For the case of the transmission to different energies each 
wave function should be normalized by the flux in the 
direction of propagation, namely, 

ttt}-* = l, (11) 

where VP is wave function and j is current operator. The 
transmission and reflection probabilities defined in such 
basis obey generalized unitarity condition which is the 
conservation of the flux during the scattering process 
from the mode with energy E to the modes E + ifku 
(see also Appendix IB] Eq. IB8|) : 

^ Re^e+Huj + Te^e+iHu = 1. (12) 

i i 

Hence we can simplify the two first terms in the expres- 
sion ((10)) : 

(13) 



L[jum] 




FIG. 7: The region of the parameter space of the length of 
the membrane and the fabrication tension (L, 7) where the 
electron motion is strongly coupled to the oscillations of the 
membrane. The region below the curve corresponds to strong 
coupling. 



We introduce the coefficients of total transmission from 
an energy E as 

^LR = E E ^E-^E+itiu' ^RL = E E ^B-IbW^) 
q i q i 

and get the expression of the generalized Landauer for- 
rnula^ 

J dE (TfME) - Ti L f R {E)) . (15) 



B. Transmission through the oscillating graphene 
membrane 

In further treatment, to get the conductivity of 
graphene with an oscillating barrier, we consider the rect- 
angular graphene sheet. First, we get the general depen- 
dence of the conductivity on the strength of the oscil- 
lations which is determined by the influence of the sub- 
strate and the applied gate voltage (the Fermi energy 
in terms of the density of the electrons on the graphene 
flake is Ep = hvF\/^n)- This dependence is not crucially 
sensitive on the shape of the graphene sheet. Second, we 
connect the conductivity of oscillating graphene to the 
distribution of the phonon frequencies (which determine 
the frequency of the oscillations in our setup) . For generic 
chaotic systems such as quantum billiards it is well known 
that the geometric shape should not influence the distri- 
bution of the eigenfrequencies after averaging over the 
ensemble (in our case, over disorder configuration). 

As we have shown in Eq. ©, the oscillations of 
the graphene membrane enter the Dirac equation which 
describes electrons in graphene as oscillating electro- 
static potential as well as pseudovector potential. Typ- 
ically, the effect of the electrostatic potential is more 
significant^, and therefore we disregard below the effect 
of the pseudomagnctic fields. The electron Hamiltonian 
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FIG. 8: Dependence of the relative change of conductance 
dG/G on the strength of the oscillations v = V e ;-ph/(fc*j) for 
fixed Fermi energy Ef (gate voltage). A: dG/G with chang- 
ing electron-phonon coupling keeping frequency of the oscil- 
lations fixed for two Fermi energies Ef/^lj — 8 (solid lines), 
Ef/FujJ = 18 (dashed lines) and three different widths W = 
10vf/u, 40vf/oj, 100vf/lo and for the length L = IOvf/co. 
Larger width corresponds to more k y harmonics taken into 
account that is why the curves with larger W are more 
smooth. B: dG/G for two Fermi energies Ef /V e i- P h = 8, 
Ef /V e i-ph = 18, fixed electron-phonon coupling, changing 
the frequency. The relative change of conductivity increases 
with increasing the strength of the oscillations. The change 
of the frequency influences on the conductivity more than the 
change of the electron-phonon coupling even though the dim- 
mensionless parameter of the strength of the oscillations is 
the same. It happens as we need to keep the values of the 
gate voltage and geometrical size in real parameters fixed. 



thus takes the form 

H = —ihvpa ■ V + Vq + V e -ph cos(wi). (16) 

This problem was considered before in a number of stud- 
ies: to investigate the direct analog of the Tien-Gordon 
effect in graphene with an oscillating lead^l, to look at 
the properties of transmission amplitudes depending on 
the angle of the incoming wave2^, and to implement 
quantum pumping via evenescent modes in graphene^. 
In our work we follow closely Ref. solving numerically 
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FIG. 9: The absolute value of the of the conductivity correc- 
tion at fixed v is approximately the same for different gate 
voltages. The large negative values of the change of the con- 
ductivity at low gate voltages are due to Klein tunneling dis- 
cussed in subsection (|IIICfl . 
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FIG. 10: Dependence of the conductivity on the gate voltage 
for a different values of the strength of the oscillations v for 
a narrow graphene strip, W = 2vf/u), L = IOvf/uj. For a 
narrow strip the conductance is quantized. Each step of the 
conductance at V s is multiplied at the values of the V 3 + hcu 
as we expect. 



the system of equations for transmission and reflection 
coefficients (as discussed in Appendix [Bj , but we go be- 
yond this work by summing up the transmission coeffi- 
cients over the different transversal momenta q to obtain 
the conductance of the system. Our final goal is to re- 
late the spectrum of phonon frequencies to the features 
of Tien-Gordon effect. 

The Tien-Gordon effect for graphene is characterized 
by the dimensionless parameter which is the ratio of the 
electron-phonon coupling V e — P h to the energy of the os- 
cillations of the phonons Huj n , 



V e -ph 

fujJn 



(17) 



Based on the expression for electron-phonon cou- 
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pling © we can estimate Vo, 

3i[eV]</> 2 [nm] 



v ~ 0.03 



y/jL[fim] 



(18) 



where we took into account the expression for the small- 
est frequency of the oscillating circular membrane Ao ~ 
2.4/L. For the case of the thermally driven oscillations 
with the amplitude of the oscillations (j> 2 ~ kbT/Muj 2 
we get 



v ~ 3 • 10" 



gi[eV]T 
7 3 / 2 L[/im] 



(19) 



with 7 characterizing a tension of the membrane, 7 = 
AL/L. For example v ~ 1 for the realistic parameters of 
graphene L = 300 11m with the tension 7 = 10~ 4 with 
the coupling potential g\ ~ 10 cV at the temperature 
T = 300 K«2S. The parameter of the strength of the oscil- 
lations can be substantially modified by the tension and 
the size of the membrane, Fig. [7] The tension depends 
on the fabrication process. While increasing the tension 
the influence of the oscillations of the membrane on the 
conductivity decreases, v decreases. With increasing the 
size of the membrane v decreases as well. Decreasing 
temperature reduces the influence of the oscillations on 
the conductivity as the coupling is proportional to tem- 
perature. We would like to note that the higher oscilla- 
tions of the membrane are less coupled to the electron 
motion, namely 



v 



,-1/2 



(20) 



We determine the influence of the oscillations of the 
potential on the conductivity using the generalized Lan- 
dauer formula (|15[) . For low temperatures, kT <C Vbi as , 
kT <C huj, we can approximate the distribution function 
in the leads by the constant, then the current through 
the system is determined by 



I(E ,v) = V bi<u T*(y) 



(21) 



The transmission coefficient is already a sum over trans- 
verse momenta and over all final energies Eq + ihuj which 
scatter from the energy Eq. Here we also took into 
account the inversion symmetry of the system setting 



rpE 

LR. 



R.L • 



For weak oscillations, v <C 1, we can make an ana- 
lytic estimation of the effect (see also Refs. I23I I24I) . The 
expression for the transmission probabilities bringing an 
electron from the energy Eq to Eq ± fru is T^ 1 (£'o) ~ 
v 2 T^(Eq)T^(E ± huj) (here and below wc use the sim- 
plified notations for the indices of the transmission co- 
efficients, namely, denote the energy by numbers which 
are i in the expression Eq + ihui). We can disregard all 
other probabilities as they are smaller in the parameter 
v. The current becomes then 



I(Eb,v)~/ (Eb)(l- 



-v 2 {Tl R {EQ+hu)+T° LR {E -?w))) , 

(22) 



where Iq is the current in the absence of the oscillations. 

For moderate values of v ~ 1 we solve the scattering 
problem numerically and then determine the conductiv- 
ity by means of the generalized Landauer formula (|15|) . 
In Fig. lIII Bl we present the dependence of the conductiv- 
ity on the strength of the electron phonon coupling v for 
different energies on the graphene flake. We notice that 
if the strength of the oscillations grows, the condictivity 
decreases and becomes independent on the width of the 
graphene sheet for large values of v. The absolute value 
of the drop of the conductivity is approximately similar 
for different values of gate voltage on graphene. 

For narrow graphene strips the conductance is quan- 
tized with increasing the gate voltage. Due to the oscilla- 
tions the steps in conductance at gate voltages V s t ep are 
reflected at the gate voltages V ste p + fiu, see Fig. [TU1 
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FIG. 11: Dependence of the conductivity on the gate voltage 
for different values of the strength of the oscillations v for a 
wide graphene strip, W = 1Qvf/w and W = AQvf/oj. The 
dependence is smooth. The conductivity slightly decreases 
(about 10 percent) at larger values of the gate voltage (See 
also Fig. [SJ. Pair creation is important at low gate volt- 
ages where conductivity significantly increases. Here the one- 
particle approximation with Dirac equation does not work and 
more precise consideration using quantum field theory should 
be applied. 
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FIG. 12: A:The area of parameters (e, v) where the Klein 
tunnelling is important. B: The region (density, size of the 
graphene sheet) where the pair creation is important obtained 
from the curve at Fig. A recalculating from the dimensionless 
units to the physics units (T = 300 K, pi = 12 eV, 7 = 10" 4 ). 



C. Klein tunneling 

The dependence of the conductivity on the gate voltage 
at high values of v (v > 1) and low gate voltages exhibits 
unexpected growth which reveals the existence of pair 
creation, see Fig. QTJ as we explain below. 

Graphene is a condensed matter material where Klein 
paradox can be measured directl y 29 ' 30 . Klein paradox 
arises in scattering of Dirac electrons on a high enough 
potential step: The tunneling probability is negative, and 
the reflection probability is greater then one. Another 
closely related counterintuitive result is that the trans- 
mission probability remains finite for arbitrary high and 
long barriers and does not depend on the length of the 
barrier. Both results come from the solution of the Dirac 
equation. It is a consequence of one-particle treatment 
of Dirac fermions. A fully consistent consideration of 
the problem requires the quantum field theory approach 
where many-particle states are taken into account. It 
was done for example in Re fs. |3lll32l and discussed in the 
context of graphene in Ref. |33|. 

Klein tunneling can be represented as tunneling with 
a creation of an electron-hole pair. When a particle-hole 
pair is created, the particle can be reflected backwards by 
the potential, which gives a positive contribution to the 



reflection probability and can make it larger than one. 
The hole moves in the same direction as the transmitted 
particle would move, thus contributing a negative term 
to the transmission probability 3 -^. In classical terms the 
motion of the hole is the same as motion of the particle 
in the opposite direction. The above picture stems from 
the single-particle (quantum-mechanical) understanding 
of the process. In fact, the particle-hole creation also 
changes the ground state of the system and introduces 
vacuum polarization. When these effects become signifi- 
cant, the calculation of the transport cannot be done in 
the simple way as above as now the potential itself in 
turn depends on the density of particles and holes. 

The equilibration of the two continua (the electron and 
the hole one) happens very fast, so that it is difficult to 
observe the pair-creation in a p-n junction directly as a 
characteristic time is very short. In a time-dependent po- 
tential the pair creation occurs as well, as it is proven for 
the general case of slowly varying time-dependent Dirac 
equation 3 -^. An attempt to consider the time-dependent 
Dirac equation in the context of graphene was done in 
Ref. HH, however, the authors did not take into account 
possible transmission to different energies. 

Obviously for low values of the gate voltage we en- 
counter Klein paradox for the time-dependent potential. 
Indeed, an incident particle with the energy E scatters 
to the energies E + nhui. For high strength of the oscil- 
lations v the scattering to the states with large n can be 
strong. It imitates the scattering on the potential step, 
with nhui playing the role of the height of the potential 
barrier. If the energy \E — nhu>\ > E then the pair cre- 
ation is allowed. The probability of the process becomes 
appreciable for sufficiently strong oscillations, v > 1. We 
identify the region of the parameters (V g ,v) where the 
conductivity achieved from the scattering problem for the 
Dirac equation shows increased conductivity due to Klein 
paradox, Fig. Q21 The increase of the transmission is a 
measurable effect, described still by one-particle Dirac 
equation for the moderate values of the tranmission in- 
crease. In the regime where transmission diverges, our 
results, based on single-particle Dirac equation, are not 
valid, and more precise quantum field theoretical treat- 
ment is needed in order to obtain correct value of the 
conductivity due to the pair-creation mechanism. Nev- 
ertheless, we can claim that the conductivity at small 
Fermi energies increases at low values of the gate voltage 
and high values of the amplitude of the oscillations. 

Thus, the dependence of conductivity of the gate volt- 
age presented at Fig. [11] is correct at large enough gate 
voltages (as it follows from Fig. IT^K). where the single 
particle approach to determine the conductivity is valid, 
the transmission is not subject to many particle process 
which involve the pair creation. We also note that the 
plots representing the dependence of the conductivity on 
the strength of the oscillations v, Fig. IIII B[ are made 
for the gate voltages large enough that the pair-creation 
process in not important. 

The increase of the conductivity due to pair creation 
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can be measured for small graphene samples where the 
coupling strength is large enough. Namely, we can ap- 
proximately connect the region on the diagram Fig. [T2l \ 
to the density of graphene n and the length of the mem- 
brane using the relation for the Fermi-energy Ep = 



G/Gq 



hvpy/rm, frequency lj '- 
strength given by Eq. (jT9 



2.4, 



— 4- and the coupling 
After these substitutions the 



region where the Klein tunneling can be significant is: 
2 • 1(T 7 0.02L[/im] v /n[10 10 cm-" 



7 3 / 2 L[/im] 
sec Fig. [Hp. 



> 0, (23) 



D. Tien-Gordon effect for several modes 

Thus far, we only considered a monochromatic excita- 
tion of the membrane. In this Subsection, we look how 
the presence of several oscillation frequencies of the mem- 
brane influences the current. If several modes are essen- 
tial, the largest contribution is given by the lowest modes, 
since for them v is the largest. The higher is the oscilla- 
tion mode the less is the correction to the conductivity is 
contributes. As we show below there is no region for pa- 
rameters in our problem where the interference between 
different frequencies is important for the conductivity, 
and the modes contribute independently: The influence 
of all modes is just the sum of individual contributions. 

In this subsection we look first at the conductivity in 
the presence of two frequencies. We solve numerically 
the scattering problem for the Hamiltonian, 

H = —ihvpa • V + Vb + V e -ph (cos(tJit) + cos(cj2^ + VO) • 

(24) 

using the same method which we set up for the case of 
one frequency. 

For weak electron-phonon coupling, vi,V2 <C 1 the 
influence of the two oscillations is additive in the lowest 
order approximation, see Eq. (|22p. Moreover, we see from 
the results of numerical treatment, displayed on Fig. [TBI 
that this additivity property is approximately valid even 
for moderate values of the oscillation strength v< ~ 1. 
Namely, the relative change of the conductivity for the 
oscillations with two frequencies (G(wi,aj 2 ) — Gq)/Gq, 
Go = G(uj = 0), equals to sum of the relative changes of 
the conductivity for the oscillations with one frequency 
(G(w) - Go)/G , 

G(wi,w 2 )-G G{oj\) — G G(w 2 ) - G 

+ ~ ■ (25) 



Go 



Go 



Go 



In the region of low energies the Klein tunneling is pre- 
served and the conductivities have the additive property, 
Fig. [TJ] The slight deviation from the additivity occurs 
since we do not take into account evanescent modes for 
determining the conductivity. 



0.98 



0.96 



0.94 



0.92 L 




- G(^,w 2 )/G 
■■■ G(l)2 — const)/Go 

■■ (G(wi) + G(w»))/G -1 



0.5 



1 



A 



1-5 _ Vd-ph 



G/G 
0.98 



0.94 



0.9 



0.86 





Ep/hwi = 8,v 2 = 1.0 _ 




















2"'"- 


G((Ji,w 2 )/Gq 


" v I , - 




v 


G(ui)/G a 


\\ 


- - - ■ (G(ui) + G(w 2 ))/G 


-i 



0.5 



1.5 



Vel-pf 
hid I 



dG/Go 



0.96 
0.94 
0.92 
0.9 




E F /hw = 8 E F jhu = 18 

v = 0.5 i 

v = 1.0 



C 



0.2 



0.4 



0.6 



0.8 



1.2 



1.4 



FIG. 13: The conductivity for the case of two modes of 
the oscillating potential. The Fermi energy is fixed. We 
assume that the modes are characterized by the same cou- 
pling to the phonons V e i- P h and different frequencies (the 
coupling strength depends on the phonon frequency, but for 
the illustration of the interaction of two frequencies this de- 
pendence is not important). Hence the dependence on v 
translates directly to the dependence on the frequency as 
1/uj. A,B,C:The conductivity dependence on frequency wi at 
fixed frequency u>2 is shown. We show G(vi, V2 = const) /Go 
and for comparison G(vi)/Go and G(v2)/Go as well as 
(G(vi) + G(v 2 ))/G - 1. We see that the effect of two fre- 
quencies is almost additive, A,B. Fig. C presents G(v) for 
different sets of parameters: at smaller energies the correc- 
tion to conductivity is more significant. 



Taking into account this property, we can determine 
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the disorder-average correction to the conductivity due 
to the presence of different phonon frequencies, 



51 



(.lis. 



\ i I dis. 

(26) 

The influence of the high-frequency modes can be eas- 
ily taken into account since their density of states is 
constant, dN/d(u> 2 ) = C, which leads to Phigh(^) = 
dN/duj = 2Cuj, and the strength of the oscillations is 
small v = V e i_ vh jhw <C 1 for high frequencies, hence the 
relative change of the conductivity for every mode is pro- 
portional to v 2 , and for the range of high energy modes 
it is 




r w."T V i tf i t ;rj " ' '""" < "" r " 



G(a)i,(o 2 )/Go 
G(a>i /G 
G(cog)/Gp 
(G(Q) 1 )+G(co 2 ))/G -T 



3 4 
E F /h co 



. (27) 

where the cut-off frequency which is the position of 
the first energy level of the oscillations, is given by the 
geometrical size of the membrane, and io max is the max- 
imum frequency of the oscillations which is determined 
by microscopic properties of the membrane at zero tem- 
perature or by kbT/H at finite temperature (we assume 
that the maximum frequency is given by comparission 
the amplitude of zero point fluctuations and the ampli- 
tude of the thermal fluctuations, namely Mu> max (/) 2 h ~ 
kbT hu rnax ). We also take into account the decrease 
of the coupling with increasing the frequency of the oscil- 
lation, Eq. (|20|) . The constant density can be estimated 
from the general consideration: The density of states of 
the oscillating membrane according to the Weyl formula 
is dN/dX 2 = (j/p )dN/duj 2 = (1 - r 2 n)L 2 /2 sa L 2 /2, 
where n and r are the number of holes and the ra- 
dius. Therefore the oscillations at high frequencies can 
decrease the conductivity by 

(^L a .,~(^)MfB- < 28 » 

We can express the change of the conductivity for the 
real sample: 

/WG\ 21 g 2 [eV]T 2 [K] / 3r[K]L[/H \ 

\ Go / dis . ~ ■ ' L 2 [ A1 m] 7 3 1, ^ J 

(29) 

The larger is the membrane the smaller is the conduc- 
tivity correction. For large values of the electron phonon 
coupling V e i-ph > one observes a noticeable decrease 
of the conductivity due to first several modes of the os- 
cillations. For example, for the graphene membrane of 
the length L = 500 nm and the tension 7 = 0.3 • 10~ 3 the 
correction is / <W ^M*) \ _ 0.015. Let me note 

\ a ° I disorder 

that the correction is very sensitive to the stretching. 



FIG. 14: Dependence of the conductivity on the Fermi energy 
in the case of a two-mode oscillating barrier. We clearly see 
that the additivity property is preserved even for low energies 
where Klein tunneling plays an important role. The deviation 
from the additivity occurs due to finite size of the system and 
lack of account of evanescent modes. 



IV. DISCUSSION 

In this Article we investigated the statistics of low- 
energy spectrum of the graphene membrane suspended 
over disordered substrate and the effects of this statistics 
on the conductivity. We modeled the graphene sheet as 
a membrane attached to randomly located circular ar- 
eas. The phonon frequencies in this system are given by 
eigenvalues of the Laplace equation with corresponding 
boundary conditions. 

It is known that the high-energy eigenvalues in this 
system obey the Wigner-Dyson statistics. This is valid 
for a generic chaotic system, where the statistics is found 
as the result of the averaging over the energy with the 
constant density of states. It was less obvious what the 
statistics of low-lying cigenfrcquencies is, since conven- 
tional analytical methods are not applicable to this situ- 
ation. For a single sample (with a given impurity config- 
uration) one can not define the statistics, since one can 
not average over energy. However, in a disordered sys- 
tem one can instead average over impurity configurations 
(positions of the attachment areas) . We performed this 
averaging numerically and found that the statistics of the 
nearest level spacing between the first five-six levels obeys 
the Wigner-Dyson probability distribution. The impor- 
tant difference between high-energy and low-energy case 
is that the density of eigenvalues is not constant for low 
energies — rather, it oscillates for the first few levels, 
which is a remnant of the existence of the positions of 
the eigenvalues of a clean circular membrane. At high 
energies, the density of eigenvalues saturates. Note also 
that the eigenvalues of a disordered membrane can never 
lie lower than the lowest eigenvalue of the clean mem- 
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brane, which creates an offset energy in the density of 
states. While changing this cut-off, in physical terms 
changing the size of the membrane, the statistics of the 
levels is not changed much. 

Our results can not be described by the non-linear a- 
model, a standard tool to deal with disordered systems, 
because in our system the disorder is not strong enough. 
We note that whereas in the non-linear cr-model effective 
bosonic fields are introduced after averaging over the dis- 
order, and these fields are confined to a sphere (compact 
sector) due to strong disorder, in our case of low-lying 
levels such fields would not be confined to the sphere. 
The indirect evidence for this is the oscillation in the 
density of levels for the first few levels. 

The phonon levels could be accessed directly, for ex- 
ample, by coupling the graphene membrane to other de- 
grees of freedom like microwave photons. However, these 
experiments presently look too involved, and therefore 
we investigated the effect of the phonons on electron 
transport through graphene. The coupling between elec- 
trons and phonons is provided by the deformation poten- 
tial. The electron transport in this situation shows the 
phonon-assistcd Tien-Gordon effect. Due to the oscilla- 
tions of the membrane with the frequency w the elec- 
tron can be transmitted from the energy E to the en- 
ergy E + ntkJ. At strong electron-phonon coupling and a 
small membrane we find that the influence of the phonon- 
assistcd Tien- Gordon effect on the conductivity can be of 
order often percent. If the frequency increases, the influ- 
ence of the phonon oscillations becomes less significant, 
and thus the lowest frequencies of the phonons have the 
most crucial influence on the change of the conductivity. 
The modes affect the current independently, and thus 
the influence of the superposition of the modes can be 
approximated by the sum of the influences of an individ- 
ual mode. 

The Tien-Gordon effect is the most pronounced in the 
case of the presence of structure in the conductivity. The 
structure, in our case the Fabry-Perot resonances, ap- 
pears in the voltage dependence of the conductivity in 
the case when tunnel barriers are located at the ends of 
the sample. Then every resonance at the energy E has its 
sattelites at the energies E + nhuj. This can be also a way 
to measure the phonon frequencies. This method can be 
applied only to the case of the narrow samples, as for 
real two-dimensional sample the Fabry-Perot resonances 
have different positions for different transverse momenta, 
and summing over the transverse momentum makes the 
Fabry-Perot peaks and their satellites indistinguishable. 

At low gate voltages and strong electron-phonon cou- 
pling the oscillation potential can create the potential 
step large enough to detect Klein tunneling process, 
namely, increase of the conductivity due to electron-hole 
pair creation. There is a range of parameters where the 
single particle Dirac equation is applicable to describe 
this process. When the reflection and transmission prob- 
abilities start to diverge, many-particle quantum field 
theoretical treatment should be performed. Unexpected 



increase of the conductivity at low gate voltages can be a 
signal of strong electron-phonon coupling. These effects 
lie outside the scope of this Article. 

Let us note that a suspended graphene membrane is 
not the only situation when low-lying eigenvalues of the 
Laplace equation can be assessed. For instance, our 
treatment also applies for a system of microwave cavities 
with randomly positioned reflective disks inside. The dis- 
tribution of photon frequencies in these cavities will have 
the same features as described in this Article. 
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Appendix A: Boundary integral method for a region 
with holes 




FIG. 15: The parameters used in Eq. (fATOl) - (|AT2l) . 



In this Appendix, we show how to solve the eigenvalue 
problem for the Laplace equation: 

A(j)(q) + Xl4>(q) = 0, for q in fi/dft (Al) 

with the boundary condition 

(j){q) = for q in d?L. (A2) 

The solution can be obtained using the boundary integral 
method widely implemented for the solution of this prob- 
lem for quantum billiards^, with the only difference that 
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we are interested in a problem with multiple boundaries, 
namely, a region with holes inside it. 

Let us first introduce the general formalism of the 
boundary integral method and then specify the problem 
for the case of several boundaries. 

We consider the normal derivative of the oscillating 
fielcU& u(s) = (d/dn)<f)(q(s)) = n(s)V0(q(s)), where s is 
a parametrization of the boundary and n is the normal 
consistent with the direction around the contour. The 
field obeys the integral equation, 



u(s) 



on 



,(s, s')u(s')ds', 



with the integral kernel 



Us,s') = -2—G n (q(s),q(st)) 
on 



(A3) 



(A4) 



where G n (q(s),q(sf)) is the Green's function of the cor- 
responding Laplace equation, 



(A + A£)G n (q,q') = *(q-q / ) . 



(A5) 



In two dimensions, the Green's function of free motion is 
Hankcl's function, Hq'K Therefore the kernel 
be expressed as 



Q n (s,s') = y cosa(s, s')H[ 1 \\ n T(s, s')), 

with t(s,s') = \q(s) — q(s')\ and 

n(q)(q( S )-q(s')) 



(A6) 



cos a(s, s') = 



t(s,s') 



For numerical implementation the equation (|A3[) is dis- 
cretized, the boundary is discretized, S(jN, and the corre- 
sponding values of uu-\ = u(s^) are taken. The condition 
that the solution exists is 



det (Q„(i,j) - 1) = 0. 



(A7) 



We vary A and monitor the left-hand side of the con- 
dition (|A7j) . thus finding the cigcnfrcqucncies of the 
Laplace equation. Let us note that numerically it is more 
efficient and more precise to look at the eigenvalues of 
Qn(i,j) ~ 1 then to look at the condition for the determi- 
nant. 

The boundary integral method can be simply general- 
ized for a region with holes. Let us represent the bound- 
ary as 90 = dflo U d£li U • • • U dil n where dtlo is outer 
boundary of the area and {f^}, i = 1,. . .,N are holes 
inside of the area. The directions of the boundaries and 
normal vectors to the holes are consistent with the direc- 
tion of the outer boundary. 

Now the function it(s) on the boundaries is represented 
as {ugn, Mani , ■ • ■ , uqq n ). For further convenience of the 
notations we use the indices O and fii instead of the pre- 
cise notations <9f2 , <9fV The system of integral equations 
is: 



->nn 



un 



ion 
Q n,n 



£ 



on i 



V n U 0.i 



(A8) 



un- 



on 



on, 



on. 



(A9) 

with Q n being the corresponding kernels given by general 
formula (lA6l). 



For example we present the kernels for our model of 
the circular outer membrane and the holes of radius r 
placed in the random positions inside of the membrane. 
We parametrize the hole position by the position of the 
center and the angle (ri,9i). The position on the outer 
membrane is given by angle 4>r, on the small circle by 4> r . 
We take the radius of the outer membrane R = 1 and for 
the holes r. All parameters are shown at Fig. [AJ By r 
we denote the distance between the points for which the 
kernel element is given, r(s, s') = |q — q'|. 



J 



• For q at the large circle and q' at small circle is (the negative sign originates from the circulation in another 
direction for the integral over the small circle) 



(s,s') = 



iXr, 1 — r cos( 



>R 



- Ti COS((j) R 



(A10) 



for q at the small circle and q' at large circle is (the negative sign originates from the direction of normal) 



i\ n r + Ti cos(0 r + 8i) 



+ 0J?) 



(Ms,*')) 



(All) 



• for both q and q' at the small circle (one negative sign originates from normal, and another one from the 
direction) 



iX n ri cos(6 t + 4> ri ) - r.j cos(4> r j + 6j) — r cos( 



(A12) 
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Appendix B: Scattering problem for the 
Tien-Gordon effect in graphene 



E + Ihu, 



We consider scattering of electrons in graphene in the 
timc-dcpcndcnt periodic potential, 



J 1 ) 

ky ^ z -k {,) ,k 



e -ikr(x e -i(E+lHu)t/ti 



H=-ihv F a-V + V + Vcos(u,t), (Bl) * t (t,s,y) = e^^fg L { J ) 



g-ifef x e -i(E+lHu)t/h 



where a = (a x ,a y ) are the Pauli matrices, Vb is the po- 
tential in the region of propagation, V is the amplitude 
and uj is the frequency of the oscillating potential. Due 
to the oscillating potential the incident wave with the en- 
ergy E scatters into waves with energies E + nhui. Our 
goal is to determine the transmission coefficients and to 
get the conductivity through the region of graphene with 
the oscillating potential. 

The time-dependent wave-function of the incident 
Dirac electron with energy E is: 



V in (t,x,y)=e ik v y 



1 



ikfjx -iEt/h 



(B2) 



with the momenta k y and k. 



(0) 



perpendicular and along 

the strip respectively, such that k 2 +(k^) 2 = E 2 , and the 
phase factor Zfc.fe = (k+ik y )/E. Reflected and transmit- 
ted wave functions are written as sums over the energies 



with the momenta k^ such as 



kl + (k { t l) ) 2 = (E + lfriu + Vb)' 2 



The wave function in the region with the time- 
dependent potential has the time dependent factor 
exp(— iV/h J cos cot) which can be expanded over the har- 
monics cxp(mwt)^, 



exp ^— i— J cosujtj = J„ exp(-inu)t) , 

(B3) 

where J„ is the Bessel function of nth order. Hence the 
wave function in the oscillating region is expressed as: 



& b (t,x,y) = e 



ik y y 



e -i(E+lHu)t/h 



A 



(0 



1 

z k il) k 



r 



B 



CO 



5> 



Hlo 



exp{—inujt) 

m 



with the momenta fe^ along the strip that obey k 2 + 

(k^) 2 = (E+Vo+lhcu) 2 . The conditions of the continuity 
of the wave function at x — 

$ in {t,0,y) + * r (t,0,y) = V b (t,0,y), 

and at x = L, 

* b (t,L,y) = V t (t,L,y), 

give the system of linear equations for the reflection and 
transmission coefficnts. 

The obtained system of equations is infinite. We solve 
the system numerically making it finite and ensure that 
the result does not depend on the number of the equation 
in the system. 

Let us note that the above system of equations takes 
into account both propagating and evanescent modes in 
the region between the leads. To determine the trans- 
mission and reflection probabilities we should take into 
account the normalization by unit flux in x direction for 
the mode of every energy, 



1 

Zk,k. 



= 1 



(B5) 



which leads us to Ne = 2k/ 'E. Note that the evanescent 
modes can not be normalized in this way. However, in 
the following we assume that the system is sufficiently 
far from the Dirac point, and the effect of evanescent 
modes on electric transport is insignificant. Then the 
transmission and the reflection probabilities are 



T, = 



k?E 



(E + lhw + V b )k 



\ti\ 2 , 



and 



Rl = 



E 



E + lhbj + V b 



In I 



(B6) 



(B7) 



The flux conservation condition (|B5[) is given by the gen- 
eralized unitarity condition (generalized as we take into 
account several modes): 



(B8) 



We ensure this property numerically to check that the 
numerical size of the system of the equations is sufficient 
to preserve the current through the system. 
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The same system of equations was considered in Rcf. 
l23l . We extend the analysis of that publication, since at 
the end we determine the full transmission probability, 
i.e. we sum over all propagating modes. We generalize 
the problem to several frequencies as well. For exam- 
ple for the case of two frequencies Ui and W2, the scat- 
tering will happen to the energies E + Ifrwi + nhu2- If 



for the treatment of the system with only one frequency 
one needs the system of the equations of the size AN, 
then for the treatment of the system with two frequen- 
cies it is AN 2 . This growth of computational complexity 
prevented us from looking at the strong coupling region 
where both parameters V/hwi t 2 are large. 
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